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Abstract 

A  fully  three-dimensional  numerical  simulation  of  plasma  density  irregularities  in  the  postsunset  equatorial 
F  region  ionosphere  leading  to  equatorial  spread  F  (ESF)  has  been  developed.  The  simulation  evolves  under 
realistic  background  conditions  including  bottomside  plasma  shear  flow  and  vertical  current.  It  has  been 
initialized  and  driven  in  accordance  with  background  measurements  obtained  from  the  Air  Force  C/NOFS 
satellite,  the  Jicamarca  Radio  Observatory,  the  Altair  radar  on  Kwajalein,  and  other  instruments.  Sim¬ 
ulations  are  seeded  with  white  noise  alone  and  are,  in  that  sense,  agnostic  with  respect  to  the  required, 
underlying  seed  mechanism.  A  combination  of  generalized  Rayleigh  Taylor  (gRT)  and  collisional  shear  in¬ 
stability  (CSI)  has  been  found  to  produce  growing  waveforms  with  key  features  that  agree  with  satellite 
and  ground-based  radar  observations  in  several  important  respects,  including  their  gross  morphology,  de¬ 
velopment  rates,  and  magnetic  signatures.  The  transient  response  of  CSI  is  consistent  with  the  observation 
of  bottomside  waves  with  wavelengths  close  to  30  km,  whereas  the  steady-state  behavior  of  the  combined 
instability  can  account  for  the  100+  km  wavelength  waves  that  predominate  in  the  F  region. 


Introduction 


The  goal  of  the  Air  Force  C/NOFS  satellite  mission  is  to  forecast  radio  propagation  disruptions  caused 
by  plasma  irregularities  formed  during  equatorial  spread  F  conditions.  The  AFRL  methodology  involves 
using  data  from  the  C/NOFS  satellite  and  elsewhere  to  initialize  large-scale  simulations  of  the  background 
equatorial  ionosphere,  let  them  evolve  in  time,  and  estimate  the  likelihood  of  instability  by  evaluating 
estimators  of  the  linear  growth  rate  of  the  generalized  Rayleigh  Taylor  (GRT)  instability.  When  and  where 
elevated  growth  rates  are  detected,  intermediate-scale  simulations  are  run  to  predict  irregularity  development, 
extent,  and  morphology  and  to  predict  radio  scintillation  levels.  Readers  may  consult  Retterer  [2005];  Retterer 
et  al.  [2005]  for  details. 

However,  sounding  rocket  experiments  on  Kwajalein,  radar  measurements  from  Altair  and  Jicamarca, 
and  numerical  simulations  from  the  Cornell  group  argue  that  an  important  intermediate  process  is  involved 
in  initiating  ESF.  At  issue  are  plasma  instabilities  that  grow  in  the  bottomside  F  and  valley  regions,  at  the 
base  of  a  stratum  of  strong  plasma  shear  where  plasma  and  neutral  drifts  are  antiparallel  (retrograde) .  The 
growth  rate  of  these  “collisional  shear”  instabilities  (CSI)  is  predicted  to  be  several  times  faster  than  GRT, 
and  the  preferred  scale  sizes  of  the  waveforms  it  produces  match  those  of  the  large-scale  waves  observed 
at  the  onset  of  ESF.  Understanding  and  modeling  these  instabilities  would  appear  to  be  a  crucial  step  in 
forecasting  spread  F. 

Under  AFOSR  support,  we  developed  a  simulation  code  that  includes  both  CSI  and  GRT  instability. 
Together,  these  processes  exhibit  rapid  growth  and  other  emergent  behavior  (characteristics  exhibited  by 
neither  one  acting  alone)  [ Aveiro  and  Hysell ,  2010].  ESF  depletions  generated  by  the  simulation  under 
realistic  background  driving  conditions  were  morphologically  similar  to  typical  Jicamarca  observations.  Sim¬ 
ulations  driven  by  C/NOFS-supplied  electric  fields  produced  ESF  plumes  that  matched,  quantitatively  and 
qualitatively,  radar  observations  from  Altair  [Aveiro  et  al.,  2011]. 

Below,  we  review  representative  C/NOFS  and  Altair  radar  observations  of  ESF  irregularities  and  show 
how  the  salient  features  of  the  observations  have  been  recovered  by  our  numerical  simulations  of  the  phe¬ 
nomena.  The  simulation  results  are  sufficiently  similar  to  the  observations  to  inspire  considerable  confidence 
regarding  their  predictive  capability. 


Representative  Radar  Observations 

The  ARPA  Long-Range  Tracking  and  Instrumentation  Radar  (ALTAIR)  is  a  VHF/UHF  radar  located  at 
Kwajalein  Atoll  (8.8°N,  167. 5°E)  [Tsunoda  et  al.,  1979].  Figure  1  displays  ALTAIR  UHF  (422  MHz)  data 
for  the  night  of  April  29,  2009.  ALTAIR  observations  made  perpendicular  to  the  geomagnetic  field  at  F 
region  heights  represent  a  combination  of  coherent  and  incoherent  scatter,  whereas  off-perpendicular  scans 
reflect  only  incoherent  scatter. 

The  top  panel  of  Figure  1  shows  ALTAIR  scans  made  pointing  perpendicular  to  the  geomagnetic  field 
around  09:45  UT  (LT  ssUT+11  hr  in  the  radar  longitudinal  sector).  By  this  time,  irregularities  had  formed 
at  the  base  of  the  F  region  as  bottom-type  layers.  These  were  patchy  and  spaced  by  about  35  km.  They 
are  common  precursors  of  fully-developed  ESF,  whereas  bottom-type  layers  tend  to  be  more  continuous 
(beam-filling)  or  absent  on  nights  when  topside  ESF  does  not  develop  [Hysell  et  al,  2004]. 

ALTAIR  scans  made  pointing  off-perpendicular  to  the  geomagnetic  field  at  10:30  UT  and  11:26  UT  are 
shown  in  the  middle  and  bottom  panels  of  Figure  1,  respectively.  By  about  10:30  UT,  the  irregularities  were 
still  confined  to  the  base  of  the  F  region,  but  the  predominant  horizontal  scale  had  increased  to  100-200 
km.  Between  about  10:30  and  11:26  UT,  the  bottomside  depletions  underwent  strong  vertical  development. 
Depleted  channels  penetrated  to  the  topside  while  retaining  their  100-200  km  horizontal  spacing.  The 
depletions  had  characteristic  tilts  and  showed  evidence  of  bifurcation  and  secondary  instabilities.  Wedges  of 
enhanced  plasma  also  appeared  to  be  driven  downward  in  between  ascending  depletions. 
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Figure  1:  ALTAIR  radar  scans  for  April  29,  2009,  showing  backscatter  power,  range  corrected  and  scaled 
to  electron  density  on  a  logarithmic  scale:  (top)  perpendicular  scan  at  09:45  UT,  off-perpendicular  scans  at 
(middle)  10:30  UT  and  (bottom)  11:26  UT.  The  rightmost  panels  depict  vertical  electron  density  profiles. 
Note  that  LT  «  UT  +  11  hr. 
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Representative  Satellite  Observations 


To  drive  simulations  of  the  April  29,  2009  ESF  event,  in  situ  electric  field  measurements  from  the  Ion  Ve¬ 
locity  Meter  (IVM)  and  Vector  Electric  Field  Instrument  (VEFI)  on  board  the  Communications/Navigation 
Outage  Forecasting  System  (C/NOFS)  satellite  were  used.  IVM  is  comprised  of  an  ion  drift  meter  and  an 
ion-retarding  potential  analyzer  that  together  provide  measurements  of  the  ion  drift  vector,  the  ion  temper¬ 
ature,  and  the  major  ion  composition.  The  VEFI  instrument  consists  of  an  electric  field  meter,  a  fluxgate 
magnetometer,  an  optical  lightning  detector,  and  a  fixed-bias  Langmuir  probe.  Among  its  products  are 
measurements  of  DC,  quasi-DC,  and  AC  electric  fields. 

Figure  2  shows  the  zonal  electric  field  (based  on  IVM  ion  drift  and  VEFI  magnetometer  measurements) 
and  (IVM)  ion  density  measurements  with  the  respective  orbital  information  tabulated  below  (geographic 
and  magnetic  coordinates,  altitude,  and  universal  and  solar  local  times).  The  results  correspond  to  C/NOFS 
passes  during  the  post-sunset  period  on  April  29,  2009,  centered  at  the  Kwajalein  longitudinal  sector.  For 
comparison,  satellite  passes  at  (top  panels)  ~  09:03  UT,  (middle  panels)  ~  10:47  UT  and  (bottom  panels) 
~  12:31  UT  are  shown. 

During  the  first  C/NOFS  pass  at  about  09:03  UT,  the  zonal  electric  field  was  essentially  zero,  and  no 
density  or  electric  field  perturbations  were  detected  in  the  altitude  range  between  495-520  km.  By  the  time 
of  the  next  pass  at  about  10:47  UT,  the  zonal  electric  field  had  increased  to  0.5-0.75  mV/m,  and  small 
electric  field  irregularities  were  emerging  between  435-450  km  altitude.  That  the  background  plasma  density 
increased  as  the  satellite  altitude  decreased  suggests  that  it  was  flying  in  the  topside,  which  is  consistent 
with  the  evidence  of  Figure  1. 

By  the  time  the  satellite  passed  through  the  Kwajalein  sector  again  at  about  12:30  UT,  the  electric 
field  irregularities  had  become  large  (up  to  about  4  mV/m).  Modest  plasma  density  irregularities  were  also 
present  at  altitudes  near  400  km.  In  view  of  Figure  1,  however,  the  satellite  most  likely  flew  just  above  the 
most  intense  plasma  density  irregularities. 


Numerical  Model  and  Simulations 

Here,  we  describe  our  numerical  scheme  for  the  simulation  of  the  plasma  density  irregularities  in  the  equatorial 
spread  F  event.  The  3-D  simulation  updates  the  plasma  density  and  electrostatic  potential  in  time  assuming 
quasineutrality  and  momentum  balance.  Inertia  is  neglected  in  the  present  incarnation.  The  characteristics 
of  the  neutral  atmosphere  are  imported  from  empirical  models.  The  background  electric  field  is  specified  and 
controls  the  forcing.  Simulations  are  performed  with  realistic  background  conditions  including  bottomside 
plasma  shear  flow  and  the  attendant  vertical  current. 

Rather  than  assuming  equipotential  magnetic  field  lines,  the  simulation  solves  the  potential  equation  in 
three  dimensions,  allowing  us  to  keep  track  of  the  distribution  of  parallel  currents  forced  by  the  dynamos 
and  of  their  closure.  This  also  accommodates  the  possible  effects  of  drift  waves  and  related  instabilities 
that  may  play  secondary  roles  in  ESF  [Huba  and  Ossakow,  1979],  [ LaBelle  et  al.,  1986],  [Hysell  et  al,  2002]. 
For  instance,  the  full  3-D  treatment  of  warm  plasma  shows  that  perturbations  on  the  surface  of  ionospheric 
plasma  clouds  propagate  and  twist  into  a  barber  pole  configuration  [Drake  and  Huba,  1987],  [Drake  et  al., 
1988],  [Zalesak  et  al,  1990].  As  a  consequence,  the  phase  of  the  corrugations  varies  strongly  along  B,  the 
plasma  cloud  loses  its  flute-like  characteristic  (i.e.,  k-B  /0),  and  the  resulting  finite  ku  has  a  dissipative 
influence  on  the  instability.  We  do  not  know  the  extent  to  which  similar  phenomenon  might  affect  the  ESF 
development,  but  a  study  is  underway.  Finally,  it  is  worth  mentioning  that  magnetic  field  lines  cease  to  be 
equipotentials  near  the  bottom  of  the  E  region,  decoupling  the  plasma  and  decreasing  the  conductive  loading 
on  the  dynamo,  which  may  influence  the  overall  circulation. 
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Figure  2:  Ion  density  (cm-3)  and  zonal  electric  field  (mV/m)  at  three  different  passes  through  the  Kwajalein 
sector:  (top  panel)  ~  0903  UT,  (middle  panel)  ~  1047  UT,  and  (bottom  panel)  ~  1230  UT.  The  tabulated 
entries  are  universal  time  (UT),  magnetic  latitude  (ML AT),  solar  local  time  (SLT),  geographic  latitude 
(GLAT),  longitude  (GLON),  and  altitude  (ALT)  in  km. 
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Conservation  Laws 


Neglecting  stress,  heat  flow,  and  inertia,  the  fluid  continuity  equation  and  momentum  equations  for  electrons 
and  ions  are  given  as: 


dna 

dt 


V  •  (n-ava)  Pa 


(1) 


0=  -  S7pa  +  naqa{  E  +  vQxB)  (2) 

a^/3 

+  mana[ g  -  van{va  -  U)  -  ^2  vap{va  -  v^)], 

P 

where  pa  =  ?zQfcsTQ  is  the  pressure.  The  variables  na,  qa,  ma,  Ta,  Pa,  La,  and  vQ  are  number  density, 
electric  charge,  mass,  temperature,  production  and  loss  rates,  and  velocity  of  the  species  a,  respectively. 
The  loss  term  includes  charge  exchange  and  dissociative  recombination  for  a  given  species.  Nighttime  pho¬ 
toionization  is  neglected,  but  production  due  to  charge  exchange  is  considered.  Here,  (3  represents  the  other 
ionized  bodies  that  exchange  momentum  or  energy  with  a.  The  terms  va$  and  van  are  the  collision  frequen¬ 
cies  with  charged  particles  and  neutrals,  respectively.  The  terms  E,  B,  U,  g,  and  ks  represent  the  electric 
field,  magnetic  field,  neutral  wind,  gravity,  and  the  Boltzmann  constant,  respectively.  Equation  1  is  not  in 
conservation  form,  since  it  includes  chemical  processes.  However,  it  can  be  converted  to  conservative  form 
through  variable  transformation.  Two  simplifications  are  applied  to  Equation  1  and  Equation  2:  (a)  the 
neglect  of  diagmagnetic  drifts  (VnQ  x  B),  since  their  associated  flux  divergence  is  very  small  when  compared 
to  the  other  terms,  and  (b)  Coulomb  collisions  (the  last  term  in  the  RHS  of  Equation  2)  are  treated  implicitly 
and  only  included  in  the  parallel  direction. 


Numerical  Scheme 


The  simulation  algorithm  performs  two  computations.  First,  the  self-consistent  electric  field  is  found  using 
a  3-D  potential  solver  based  on  quasineutrality,  which  can  be  expressed  as 


E  •  V<F 


=  V 


E  •  (E0  +  U  x  B)  +  D  ■  V?r  - 


(3) 


where  E0  is  the  background  electric  field,  E  is  the  conductivity  tensor,  D  includes  the  diffusivity  tensor, 
5  ■  g  is  gravity-driven  current  density,  and  d>  is  the  electrostatic  potential  that  must  arise  to  preserve  a 
solenoidal  current  density.  Explicit  definitions  of  the  tensors  are  given  in  Shume  et  al.  [2005].  Note  that 
the  resulting  electric  field  is  the  combination  of  the  background  and  electrostatic  fields  (i.e. ,  E  =  E0  —  Vd>). 
Equation  3  involves  a  partial  differential  equation  for  the  potential  in  three  dimensions.  It  is  solved  using 
the  BiConjugate  Gradient  Stabilized  (BiCGStab)  method  (e.g.  van  der  Vorst  [1992])  using  the  algorithms 
described  by  Saad  [1990]. 

Second,  we  solve  a  discretized  version  of  Equation  1  for  each  ion  species  using  a  monotone  upwind 
scheme  for  conservation  laws  (MUSCL)  (a  pedagogical  review  can  be  found  in  Trac  and  Pen  [2003])  directly 
applicable  to  the  ion  continuity  problem  (once  the  recombination  term  has  been  absorbed  into  a  redefinition 
of  na  so  that  Equation  1  becomes  a  conservation  equation,  as  necessary).  The  characteristic  of  the  neutral 
atmosphere  (densities,  temperature,  and  wind  velocity)  are  updated  in  time  on  the  basis  of  inputs  from 
climatological  models.  The  background  electric  field  is  specified  and  partly  controls  the  forcing.  The  ion 
composition  includes  0+,  NO+,  and  0^".  The  model  was  constructed  using  tilted  magnetic  dipole  coordinates 
( p ,  q ,  <j> ),  where  the  tilt  is  matched  to  the  magnetic  declination  in  the  longitude  of  interest.  In  our  terminology, 
p  represents  the  Mcllwain  parameter  (L),  q  is  the  magnetic  co-latitude,  and  </>  is  longitude  (see  e.g.,  Hysell 
et  al.  [2004]). 
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Empirical  Model  Drivers 

To  initialize  the  model  runs,  we  derive  plasma  number  densities  from  the  Parametrized  Ionospheric  Model 
(PIM),  a  parametrization  of  the  output  of  several  regional  theoretical  model  outputs  generated  for  different 
climatological  conditions  and  tuned  somewhat  to  agree  with  data  from  Jicamarca  and  elsewhere  [Daniell 
et  al.,  1995].  PIM  reproduces  electron  density  profiles  from  Kwajalein  and  Jicamarca  fairly  closely  if  scaled 
slightly  to  account  for  day-to-day  variability  [Aveiro  and  Hysell ,  2010].  The  background  number  density 
profiles  used  for  initial  conditions  are  intended  to  reproduce  measurements  from  ALTAIR,  including  the 
F  region  peak  height  and  peak  density,  the  bottomside  density  scale  height,  and  the  density  in  the  valley 
region. 

Plasma  mobilities  are  calculated  using  neutral  composition  and  temperature  estimates  from  the  Mass 
Spectrometer  and  Incoherent  Scatter  (NRL-MSISE00)  model  [Picone  et  al.,  2002]  and  ionospheric  composi¬ 
tion  estimates  from  the  International  Reference  Ionosphere  (IRI-2007)  model  [Bilitza  and  Reinisch,  2008].  We 
take  the  ionosphere  and  neutral  atmosphere  to  be  in  thermodynamic  equilibrium  after  sunset  ( Tn  =  Te  =  Ti). 
Expressions  for  the  ion-neutral  and  electron-neutral  collision  frequencies  used  to  compute  conductivities  can 
be  taken  from  Richmond  [1972].  Expressions  for  the  Pedersen,  Hall,  and  direct  mobilities  and  diffusivities 
themselves  are  found  in  Kelley  [2009],  for  example. 

For  our  runs,  the  zonal  neutral  winds  are  obtained  from  the  Horizontal  Wind  Model  (HWM-07)  [Drob, 
2008] .  As  with  PIM,  these  may  be  tuned  slightly  to  account  for  day-to-day  variability.  The  model  is  sensitive 
to  the  wind  speeds  throughout  the  thermosphere,  since  the  growth  rate  of  the  collisional  shear  instability 
depends  on  the  counter  streaming  flow  of  the  neutrals  and  plasma.  Meridional  and  vertical  winds  were 
neglected  in  the  simulations  presented  here. 

Simulation  Initialization 

Our  simulation  was  cast  on  a  rectangular  grid  139x133x189  points  wide  in  ( p,q,(f> )  space  constructed  using 
a  tilted  magnetic  dipole  coordinate  system.  A  cut  through  the  equatorial  plane  spans  altitudes  between 
90-510  km  and  longitudes  between  ±  6°.  The  flux  tubes  covered  by  the  parallel  coordinate  all  reach  to 
the  lower  E  region.  The  runs  shown  in  this  section  are  centered  on  the  dip  equator  near  Kwajalein  (5.5°  N 
latitude,  166.5°  E  longitude,  7.4°  declination)  with  background  conditions  modeled  for  April  29,  2009.  To 
seed  the  simulation  run,  we  added  independent  Gaussian  white  noise  to  the  initial  number  density  with  a 
20%  relative  amplitude.  The  local  time  dependence  of  the  background  zonal  electric  field  ( E was  estimated 
using  measurements  from  IVM  and  VEFI  on  board  the  C/NOFS  satellite.  E $  was  taken  to  be  constant  in 
the  F  region  but  its  amplitude  is  reduced  below  the  valley  to  avoid  the  growth  of  electrojet  gradient-drift 
instabilities.  Based  on  those  measurements  E $  was  modeled  as  0  mV/m  before  21:00  LT,  and  0.75  mV/m 
after  22:00  LT  with  a  gradual  transition  using  a  hyperbolic  tangent  function  in  between.  For  example,  £^,=34 
/zV/m  on  the  western  boundary  and  0.75  mV/m  on  the  eastern  boundary  of  the  simulation  by  10:40  UT, 
since  those  locations  represent  21:22  LT  and  22:10  LT,  respectively. 


Simulation  Results 

The  initial  conditions  of  our  simulation  at  to=9:30  UT  («20:30  LT  for  the  Kwajalein  sector)  are  shown 
in  Figure  3.  The  top  panel  shows  the  background  plasma  density.  The  initial  conditions  are  consistent 
with  ALTAIR  radar  observations  for  April  29,  2009,  at  that  time.  The  middle  panel  of  Figure  3  shows  the 
perpendicular  current  density  with  the  magnitude  and  direction  indicated  by  the  color  disc  in  the  lower 
panel.  There  are  two  major  components  to  the  current  density  that  are  combined  in  this  depiction.  First,  a 
zonally  eastward  current  driven  by  gravity  exists  throughout  the  region,  showing  a  maximum  at  the  F  peak 
near  330  km  altitude.  In  addition,  a  zonally  eastward  F  region  wind  drives  a  vertical  upward  current  in  the 
region  above  about  200  km.  The  empirically  modeled  wind  is  almost  constant  with  altitude,  and  thus  the 
vertical  current  scales  like  the  background  ion  density,  showing  a  peak  near  330  km  and  a  strong  vertical 
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gradient  in  the  bottomside  between  280  km  and  240  km  altitude.  Below  150  km,  the  current  is  influenced 
by  both  the  local  zonal  winds  and  the  enhancement  of  the  zonal  current  in  the  electrojet. 

The  divergence  in  the  vertical  current  density  is  supported  by  a  vertical  electric  field  shown  to  the  right 
of  the  middle  panel.  Above  the  bottomside,  the  field  is  directed  downward,  consistent  with  an  eastward 
E  X  B  drift  in  the  same  direction  as  the  dynamo  wind.  Across  the  bottomside,  the  vertical  electric  field 
changes  direction,  producing  a  sheared  zonal  E  x  B  drift,  which  provides  free  energy  for  the  CSI  below  the 
shear  node  due  to  the  differential  velocity  between  plasma  and  neutrals. 

The  bottom  panel  of  Figure  3  shows  equatorward  currents  in  both  hemispheres  on  field  lines  that  thread 
the  bottomside  in  the  altitude  range  from  250  km  to  280  km.  These  currents  feed  the  upward  current  driven 
by  the  F  region  zonal  wind  seen  in  the  middle  panel.  Above  the  bottomside,  the  meridional  current  is 
distributed  across  all  the  field  lines  and  is  poleward  in  both  hemispheres,  serving  to  close  the  current  loop 
driven  by  the  F  region  dynamo. 

Figure  4  shows  the  results  for  times  matching,  roughly,  the  three  ALTAIR  scans  and  C/NOFS  passes. 
Note  that  the  vertical  axes  indicate  the  apex  altitude,  i.e.,  the  altitude  mapped  to  the  magnetic  equator. 
For  example,  the  magnetic  field  line  that  crosses  over  Kwajalein  at  380  km  altitude  maps  to  400  km  altitude 
at  the  magnetic  equator. 

The  top  panel  of  Figure  4  depicts  the  evolution  of  ionospheric  irregularities  in  the  absence  of  a  background 
electric  field  due  entirely  to  the  collisional  shear  instability.  The  spatial  scales  of  the  irregularities  forming  at 
the  base  of  the  bottomside  range  between  30-50  km.  These  irregularities  cannot  undergo  vertical  development 
and  remain  confined  to  altitudes  where  the  plasma  flow  is  retrograde. 

By  the  time  of  the  second  panel  in  Figure  4,  the  transient  phase  of  the  CSI  was  ending,  and  the  asymptotic 
phase  was  underway.  The  30-50  km  waves  at  the  base  of  the  bottomside  had  coalesced  into  large-scale  waves 
with  scale  sizes  between  100-200  km.  Under  the  action  of  an  increasing  background  zonal  electric  field,  the 
bottomside  irregularities  begin  to  exhibit  vertical  development,  more  so  on  the  eastern  side  of  the  simulation 
volume,  which  had  experienced  an  eastward  background  electric  field  for  more  time. 

In  the  final  panel  of  the  simulation,  under  the  influence  of  a  significant  background  electric  field,  the 
irregularities  penetrate  to  the  topside  where  gravity-driven  currents  can  contribute  further  to  growth.  Some 
of  the  original  intermediate  scale  structure  from  the  first  panel  in  the  simulation  survive,  contributing 
to  the  overall  fine  structure.  Bifurcation  is  evident,  as  well  as  secondary  instabilities  growing  mainly  on 
the  westward  walls  of  the  primaries.  Some  descending,  wedge-shaped  plasma  enhancements  accompany  the 
ascending  depletions.  The  characteristic  westward  tilt  of  the  irregularities  is  associated  with  the  combination 
of  vertical  and  zonal  current  and  the  combined  CSI  and  gRT  instabilities.  As  the  simulation  was  initialized 
with  white  noise,  and  as  gRT  is  a  universal  instability  (nearly  scale-size  independent),  the  scale  sizes  of  the 
irregularities  evident  here  are  a  consequence  of  CSI. 

Had  a  satellite  passed  through  the  simulated  ionosphere  at  379  km  altitude  over  the  Kwajalein  sector 
(~400  km  apex  altitude)  at  about  11:13  UT,  the  electron  density  and  electric  field  perturbations  would 
have  appeared  as  shown  in  Figure  5.  Plasma  depletions  are  present  preferentially  in  the  eastern  half  of  the 
simulation  grid.  The  western  walls  of  the  irregularities  have  steeper  density  gradients  than  the  eastern  walls. 
Perturbations  in  the  zonal  electric  fields  of  the  order  3-5  mV/m  are  present,  collocated  with  large  and  small 
plasma  depletions  alike.  Since  the  pass  is  a  few  kilometers  above  the  top  of  most  of  the  plasma  depletions, 
relative  perturbations  in  the  density  are  small  in  all  but  one  case,  but  the  fringing  fields  penetrate  to  higher 
altitudes  and  are  detectable  in  all  cases. 


Overall  Findings 

In  simulation,  the  combination  of  CSI  and  gRT  produced  growing  waveforms  with  characteristics  that 
matched  the  observations  made  by  ALTAIR  radar.  The  dominant  plasma  irregularities  with  (initially) 
decakilometric-scale  took  approximately  lOOmin  to  evolve  from  the  bottomside  up  to  the  F  peak.  By  the 
time  they  reached  the  F  peak,  they  had  coalesced  into  waveforms  with  100  km+  scale  sizes.  This  behavior 
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Figure  3:  Initial  conditions  for  ESF  model  run.  Top:  Plasma  number  density.  Middle:  Transverse  current 
density  in  the  equatorial  plane  with  equipotential  contours  (in  kV)  superimposed  in  white.  Middle  right: 
vertical  electric  field.  Bottom:  current  density  in  the  meridional  plane.  The  current  densities  are  vector 
quantities  with  magnitudes  and  directions  indicated  by  the  color  wheel  at  the  bottom  right.  The  maximum 
scale  is  3  nA/m2  in  the  equatorial  plane  and  30  nA/m2  in  the  meridional  plane.  Note  that  diamagnetic 
currents  are  not  represented  in  this  diagnostic  figure. 
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Figure  4:  Simulated  plasma  densities  in  a  cut  through  the  Kwajalein  latitudinal  plane.  The  three  panels 
show  results  for  simulation  times  of  0h35m,  lhlOm,  and  lh43m,  respectively,  after  a  start  time  t0  =  9:30  UT 
(LT  =  UT  +  11  hr  at  the  horizontal  center  of  the  simulation). 
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Figure  5:  Electron  density  (cm  3)  and  zonal  electric  field  (mV/m)  at  11:13  UT  obtained  from  a  longitudinal 
cut  through  the  simulation  at  379  km  altitude  over  the  Kwajalein  sector  (~400  km  apex  altitude). 
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was  observed  in  ALTAIR  scans  and  simulation  results. 

The  transient  response  of  CSI  explained  ~  30-50  km  wave  observations,  a  feature  commonly  observed  by 
backscatter  scatter  radars  in  the  form  of  patchy  bottom-type  ESF.  These  intermediate-scale  irregularities 
were  observed  initially  in  both  the  simulations  and  the  radar  observations  but  did  not  reappear  after  they 
coalesced  into  emergent  100  km+  scale  irregularities.  This  result  agrees  with  observations  made  at  other 
longitudinal  sectors,  e.g.,  at  Jicamarca  [ Hysell ,  2000]. 

During  the  asymptotic  growth  phase,  when  intermediate-scale  depletions  dominated  the  scenario,  some  of 
the  irregularities  remained  confined  to  the  bottomside,  while  others  experienced  more  vertical  development. 
Had  these  irregularities  been  observed  with  coherent  scatter,  e.g.  by  the  Jicamarca  radar,  they  would  have 
been  labeled  bottomside  ESF.  The  irregularities  in  the  middle  panels  of  the  simulation  (Figure  4b)  and  radar 
data  (Figure  lb)  are  comparable  in  terms  of  scale  sizes  and  altitude  spans. 

The  last  stage  of  development  of  ESF  occurred  when  the  vertical  extension  of  the  plasma  depletion  was 
large  enough  that  the  irregularities  crossed  through  the  peak  into  the  topside.  Similar  features  are  observed 
in  both  the  simulation  (Figure  4c)  and  the  radar  data  (Figure  lc).  Depletions  that  reached  the  F  peak  were 
spaced  by  hundreds  of  km  and  displayed  tilts  close  to  30°  westward  from  zenith. 

Measurements  made  by  IVM  showed  that  electric  field  perturbations  of  the  order  of  a  few  mV/m  were 
sometimes  observed  even  when  depletions  in  the  densities  were  not  detected.  Similar  features  were  present 
in  the  simulation  and  seemed  to  be  a  natural  characteristic  of  the  ESF  morphology.  The  simultaneous 
measurements  of  small  relative  perturbations  in  the  density  and  large  zonal  electric  field  perturbations 
seems  to  be  a  consequence  of  a  satellite  pass  near  the  top  of  a  plasma  depletion.  Since  the  depletion  itself 
was  not  intercepted,  it  was  not  observed,  but  the  fringing  fields  penetrated  to  higher  altitudes  and  were 
detected. 

Zonal  neutral  winds  in  the  bottomside  and  valley  regions  exerted  control  over  initial  irregularity  formation. 
This  sensitivity  to  the  winds  was  expected,  since  the  growth  rate  of  the  collisional  shear  instability  depends 
on  the  counter  streaming  flow  of  the  neutrals  and  plasma  [ Hysell  et  al. ,  2006].  As  wind  measurements  in  the 
lower  thermosphere  and  mesosphere  show  considerable  variability  [Larsen,  2002],  sensitivity  tests  to  winds 
at  those  altitudes  need  to  be  performed  to  evaluate  their  effects  on  ESF  initiation. 


Conclusions 

The  combination  of  generalized  Rayleigh  Taylor  and  collisional  shear  instability  produced  growing  waveforms 
with  characteristics  matching  satellite  and  radar  observations  in  terms  of  their  gross  morphology  and  rates  of 
development.  Growing  intermediate-  (~30  km)  and  large-scale  (100+  km)  waves  resulted  from  the  nonlocal 
and  non-normal  characteristics  of  the  system  and  were  not  imposed  a  priori.  The  transient  response  of  CSI 
explained  rapid  intermediate-scale  wave  emergence  at  the  base  of  the  F  region,  whereas  the  steady-state 
behavior  seemed  to  account  for  the  large-scale  waves  that  ultimately  dominated.  The  three  different  types  of 
ESF  (bottom-type,  bottomside,  and  topside)  appeared  to  be  a  result  of  the  time  sequencing  of  the  transient 
and  steady-state  response  of  CSI  in  relation  to  the  zonal  electric  field. 

A  major  objective  of  this  research  was  to  simulate  the  electrodynamics  of  ionospheric  plasma  irregularities, 
including  those  in  ESF,  in  the  direction  parallel  to  the  geomagnetic  field.  The  results  of  ESF  simulation 
presented  here  suggests  that  a  full  3-D  treatment  is  necessary  in  order  to  track  the  distribution  of  parallel 
currents  and  the  parallel  closure  of  currents  forced  by  dynamos.  Parallel  currents  are  found  to  have  effects 
on  the  overall  circulation  and  contribute  to  ESF  morphology.  Studies  are  underway  to  quantify  the  effects 
and  mechanisms  at  work.  Parallel  currents  are  also  necessary  to  account  for  spaceborne  magnetometer 
measurements,  such  as  those  produced  routinely  by  the  CHAMP  satellite  during  ESF  encounters.  Such 
measurements  complement  the  kind  of  in  situ  measurements  provided  by  C/NOFS  and  provide  a  rather 
complete  specification  of  the  irregularities  and  associated  electrodynamics  that  constitute  ESF.  Our  fully 
three-dimensional  numerical  model  and  analysis  of  ESF  appears  to  be  able  to  reproduce  all  of  the  salient 
features  in  the  combined  satellite  and  ground-based  radar  datasets. 
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